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Abstract 

A  framework  for  solving  variational  problems  and  partial  differential  equations  that  define  maps 
onto  a  given  generic  manifold  is  introduced  in  this  paper.  We  discuss  the  framework  for  arbitrary  target 
manifolds,  while  the  domain  manifold  problem  was  addressed  in  [3].  The  key  idea  is  to  implicitly 
represent  the  target  manifold  as  the  level-set  of  a  higher  dimensional  function,  and  then  implement  the 
equations  in  the  Cartesian  coordinate  system  of  this  new  embedding  function.  In  the  case  of  variational 
problem,  we  restrict  the  search  of  the  minimizing  map  to  the  class  of  maps  whose  target  is  the  level- 
set  of  interest.  In  the  case  of  partial  differential  equations,  we  implicitly  represent  all  the  equation 
characteristics.  We  then  obtain  a  set  of  equations  that  while  defined  on  the  whole  Euclidean  space,  they 
are  intrinsic  to  the  implicit  target  manifold  and  map  into  it.  This  permits  the  use  of  classical  numerical 
techniques  in  Cartesian  grids,  regardless  of  the  geometry  of  the  target  manifold.  The  extension  to  open 
surfaces  and  submanifolds  is  addressed  in  this  paper  as  well.  In  the  latter  case,  the  submanifold  is  defined 
as  the  intersection  of  two  higher  dimensional  surfaces,  and  all  the  computations  are  restricted  to  this 
intersection.  Examples  of  the  applications  of  the  framework  here  described  include  harmonic  maps  in 
liquid  crystals,  where  the  target  manifold  is  an  hypersphere;  probability  maps,  where  the  target  manifold 
is  an  hyperplane;  chroma  enhancement;  texture  mapping;  and  general  geometric  mapping  between  high 
dimensional  surfaces. 
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1  Introduction 


In  a  number  of  applications  in  mathematical  physics,  image  processing,  computer  graphics,  and  medical 
imaging,  we  have  to  solve  variational  problems  and  partial  differential  equations  defined  on  a  general  man¬ 
ifold  Af  {domain  manifold),  mapping  the  data  onto  another  general  manifold  J\f  {target  manifold).  That  is, 
we  deal  with  maps  from  M.  to  J\f.  When  these  manifolds  are  for  example  three  dimensional  surfaces,  the 
implementation  of  the  corresponding  gradient  descent  flow  or  the  given  PDE’s  is  considerably  elaborated. 
In  [3]  we  have  shown  how  to  address  this  problem  for  general  domain  manifolds,  while  restricting  the  target 
manifolds  N  to  the  trivial  cases  of  the  Euclidean  space  or  hyper-spheres.  The  key  idea  was  to  implicitly 
represent  the  domain  surface  as  the  (zero)  level-set  of  a  higher  dimensional  function  f,  and  solve  the  PDE 
in  the  Cartesian  coordinate  system  of  this  new  embedding  function.  The  technique  was  justified  and  demon¬ 
strated  in  [3].  It  is  the  goal  of  this  paper  to  show  how  to  work  with  general  target  manifolds,  and  not  just 
hyper-planes  or  hyper-spheres  as  previously  reported  in  the  literature.  Inspired  by  [3],  we  also  embed  the 
target  manifold  J\f  as  the  (zero)  level-set  of  a  higher  dimensional  function  fj.  That  is,  when  solving  the 
gradient  descent  flow  (or  in  general,  the  PDE),  we  guarantee  that  the  map  receives  its  values  on  the  zero 
level-set  of  'f.  The  map  is  defined  on  the  whole  space,  although  it  never  receives  values  outside  of  this 
level-set.  Examples  of  applications  of  this  framework  include  harmonic  maps  in  liquid  crystals  {J\f  is  an 
hypersphere)  and  3D  surface  warping  [43].  In  this  last  case,  the  basic  idea  is  to  find  a  smooth  map  between 
two  given  three  dimensional  surfaces.  Due  to  the  lack  of  the  new  frameworks  introduced  here  and  in  [3], 
this  problem  is  generally  addressed  in  the  literature  after  an  intermediate  mapping  of  the  surfaces  onto  the 
plane  is  performed  (see  also  [25,  46]).  With  these  novel  frameworks,  direct  three  dimensional  maps  can 
be  computed  without  any  intermediate  mapping,  thereby  eliminating  their  corresponding  geometric  distor¬ 
tions  [31].  Eor  this  application,  as  in  [43],  boundary  conditions  are  needed,  and  how  to  add  them  to  the 
frameworks  introduced  here  and  in  [3]  is  addressed  in  [31]. 

To  introduce  the  ideas,  in  this  paper  we  concentrate  on  flat  domain  manifolds.  ^  When  combining  this 
framework  with  the  results  on  [3],  we  can  of  course  work  with  general  domains  and  then  completely  avoid 
other  popular  surface  representations,  like  triangulated  surfaces.  We  are  then  able  to  work  with  intrinsic 
equations,  in  Euclidean  space  and  with  classical  numerics  on  Cartesian  grids,  regardless  of  the  geometry 
of  the  involved  domain  and  target  manifolds.  In  addition  to  presenting  the  general  theory,  we  also  address 
the  problem  of  target  submanifolds  and  open  surfaces.  A  number  of  theoretical  results  complement  the 
algorithmic  framework  here  described. 

The  implicit  representation  of  surfaces  here  introduced  for  solving  variational  problems  and  PDE’s  is 
inspired  in  part  by  the  level-set  work  of  Osher  and  Sethian  [33].  This  work,  and  those  that  followed  it, 
showed  the  importance  of  embedding  deforming  surfaces  in  higher  dimensional  functions,  obtaining  more 
robust  and  accurate  numerical  algorithms  (and  topological  freedom).  Note  that  in  contrast  with  the  level-set 
approach  of  Osher  and  Sethian,  our  target  manifold  is  fixed,  what  is  “deforming”  is  the  dataset  being  mapped 
onto  it. 

Numerical  schemes  that  solve  gradient  descent  flows  and  PDE’s  onto  generic  target  manifolds  Af  (and 
spheres  or  surfaces  in  particular)  will  in  general  move  the  points  outside  of  J\f  due  to  numerical  errors.  The 
points  will  then  need  to  be  projected  back,^  see  for  example  [1,9]  for  the  case  of  Af  being  a  sphere  (where 
the  projection  is  trivial,  just  a  normalization).  Eor  general  target  manifolds,  this  projection  means  that  for 
every  point  p  E  {Af  C  iR^*)  we  need  to  know  the  closest  point  to  p  in  Af.  This  means  knowing  the 

^For  completeness,  we  will  present  the  general  equations  for  both  generic  domain  and  target  manifolds  at  the  end  of  the  paper. 
These  equations  are  easily  derived  from  [3]  and  the  work  presented  in  this  paper. 

^For  particular  flat  target  manifolds  as  the  whole  space  or  as  those  in  [34],  the  projection  is  not  needed.  Other  authors,  e.g., 
[6,  26],  have  avoided  the  projection  step  for  particular  cases,  while  in  [48]  the  authors  modify  the  given  variational  formulation  to 
include  the  projection  step. 
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distance  from  every  point  p  €  IR^  to  M  (or  at  least  all  points  in  a  band  of  J\f).  This  is  nothing  else  than  an 
implicit  representation  of  the  target  M,  being  the  particular  embedding  a  distance  function.  This  presents  an 
additional  justification  for  the  framework  here  introduced. 

In  a  number  of  applications,  the  surfaces  are  already  given  in  implicit  form,  e.g.,  [5],  therefore,  the 
framework  introduced  in  this  paper  it  is  not  only  simple  and  robust,  but  it  is  also  natural  in  those  applications. 
On  the  other  hand,  not  all  surfaces  (manifolds)  are  originally  represented  in  implicit  form.  When  the  target 
manifold  J\f  is  simple,  like  hyper-spheres  in  the  case  of  liquid  crystals,  the  implicitation  process  is  trivial. 
For  generic  surfaces,  we  need  to  apply  an  algorithm  that  transforms  the  given  explicit  representation  into 
an  implicit  one.  Although  this  is  still  a  very  active  area  of  research,  many  very  good  algorithms  have  been 
developed,  e.g.,  [14,  18,  27,  45]. 

2  The  Framework 

From  now  we  assume  that  the  target  manifold  M  is  given  as  the  zero  level  set  of  a  higher  dimensional 
embedding  :  IR^  M,  which  we  consider  to  be  a  distance  function  (this  mainly  simplifies  fhe  nofafion). 
For  fhe  case  where  AC  is  a  surface  in  fhree  dimensional  space  for  example,  fhen  ip  :  ]R^  IR.  We 
also  assume  fhaf  fhe  domain  manifold  Ad  is  llal  and  open  (as  menfioned  in  fhe  infroducfion,  general  domain 
manifolds  were  addressed  in  [3]).  We  illusfrale  fhe  basic  ideas  wifh  a  functional  from  fhe  fheory  of  harmonic 
maps.  This  is  Jusf  a  particular  example  (and  a  very  imporfanf  one),  and  from  if  will  be  clear  how  fhe  same 
argumenfs  can  be  applied  fo  any  given  variafional  problem  and  PDF.  In  parficular,  if  can  be  applied  fo 
common  Navier-Sfokes  flows  used  in  brain  warping  [31]. 

2.1  Variational  Formulation 

We  search  for  necessary  conditions  for  fhe  funclional  E[u],  defined  by 

E[u\  =  [  e[u]  dMV  (1) 

Jm 

where 

ep]  =  ^llJalll-  (2) 

fo  achieve  a  minimum.  Here,  ||  •  ||^  =  Yliji')ij  •^he  norm  of  Frobenius  and  is  fhe  Jacobian  of  fhe 
map  u  :  Ai  ^  {ip  =  0}.  Nofe  fhaf  here  we  are  already  resfricfing  fhe  map  fo  be  onfo  fhe  zero  level-sef 
of  Ip,  fhaf  is,  onfo  fhe  surface  of  inferesf  AC  (fhe  largel  manifold).  This  is  whaf  permifs  us  fo  work  wifh  fhe 
embedding  funclion  and  fhe  whole  space,  while  guaranfeeing  fhaf  fhe  map  will  always  be  onto  fhe  largel 
manifold,  as  desired.  We  use  to  nofe  fhaf  for  fhe  mosl  general  case,  fhe  funclion  is  vecforial.  Once  again, 
Ihis  energy  will  be  used  Ihroughouf  Ihis  paper  to  exemplify  our  framework.  If  will  be  clear  after  developing 
Ihis  example  fhaf  fhe  same  argumenfs  work  for  olher  variafional  formulations,  as  well  as  for  generic  PDF’s 
defined  onfo  generic  surfaces. 

Assume  fhaf  ft  is  a  map  minimizing  E{-).  Given  f  >  0,  we  conslrucl  fhe  variation 

^  A  ^ 
vt  =  u  +  t  r 

where  r  is  a  compacl  map  C°°  in  Ai.  For  an  arbilrary  x  €  A4,  we  will  in  general  nol  oblain  fhaf  Vt(x)  € 
{ip  =  0},  fhaf  is,  ip{vt{x))  0.  Therefore,  Ihis  variation  is  nol  admissible.  On  fhe  olher  hand,  we  can  from 
if  conslrucl  an  admissible  variation  via 
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Wt  = 

where  ■  R‘^  ^  O’  =  0}  is  the  projection  operator  onto  =  0}.  Note  that  sinee  ■)/)  is  a  signed 

distanee  funetion,  we  ean  simply  write  this  projeetion  operator  onto  =  0}  as 


n{^=o}(«)  =  a  -  V’(«)  VV’(a). 


Let’s  now  define 


A 


£{t)  =  E[wt\ 

Sinee  the  energy  aehieves  a  minimum  for  t  =  0, 

A  dE{t) 


fn  = 


dt 


=  0. 


(t=o) 


Let’s  eompute  this  first  variation.  We  have  that 


fo  = 


■  a  f  9yji\ 
dwl^[dxj) 

dxj  dt 


dMV 


t=o 


Moreover  (H^  stands  for  the  Hessian  of 
dwt 


dx-i 


.  f  du  dr 

-  +  i— 


and  we  observe  that 


dwf 


dxj 


(t=o) 


-  fvV'(5)  •  ^  )  VV'(5) 


dx-. 


sinee  'tp{u)  =  0.  We  ean  further  simplify  this  observing  that  0  =  =  V'tp{u)-  Therefore, 


dwf 


dxj 


_  du 
(t=0)  ^^3 


(3) 


(4) 


(5) 


With  a  bit  of  further  simple  analysis  we  ean  eompute  the  additional  derivative, 


x(£f) 

dt  dxj 


This  ehange  in  the  order  of  derivatives  is  done  in  order  to  immediately  evaluate  the  result  at  f  =  0,  thereby 
simplifying  the  following  derivative.  Following  in  an  similar  form,  we  obtain 


dwl 

dt 


=  r-  {V‘tl){wt)  ■  r)  V't/;(wt)  -  'tl’im)  H^(tyt)r 


and 


dwj 

dt 


(t=o) 


=  r  —  {'V'tp{u)  •  r)  V'tp{u). 


(6) 


(V) 
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Combining  the  above  eomputations  all  together  we  obtain 


Following  from  (3)  we  have  that^ 


(8) 


—  (r  •  V'ip{u)) 


^0  -  X] 


f 

1  dxj  dt 


^  f  (  dr  du  r 


Now,  applying  the  divergenee  theorem  we  eonelude  the  eomputation.  We  first  write 


Vr*  •  Vrr* 


and  then  apply  the  faet  Vr*  •  Vrt*  =  V  •  (r*Vtt*)  —  r*  Art*,  together  with  the  divergenee  theorem,  to  obtain 
(n  stands  for  the  outward  unit  normal  to  dM). 

X]  [  ^  X]  /  r^^dS  -  f  r^Au'dMV  (10) 

^  Jm  i  JdM  Jm 

tj  t 

To  eonelude  we  put  together  this  last  expression  with  (8),  and  after  some  algebra  we  obtain  that  £q  is 
equal  to 


,ndS-lj.[ 


AB+  j  VV'WidA,. 


The  boundary  eondition  is  eliminated  sinee  the  support  of  r  is  eompaetly  ineluded  in  Af.  To  eliminate 
the  additional  term  for  an  arbitrary  r  we  must  impose 


‘^“+(EH*[g.g])v^(«)=0.  (.2) 

This  gives  the  eorresponding  Euler-Lagrange  for  the  given  variational  problem.  Note  onee  again  from 
our  eomputations  that  in  spite  that  all  the  terms  “live”  in  the  Euelidean  spaee  embedding  the  target  manifold, 
u  will  always  map  onto  the  level-set  of  interest,  {-ip  =  0},  and  therefore,  onto  the  surfaee  of  interest.  This 
is  guaranteed  by  this  equation,  no  additional  eomputations  are  needed.  This  is  the  beauty  of  the  approaeh, 
while  working  freely  on  the  Euelidean  spaee  (and  therefore  with  Cartesian  numeries),  we  ean  guarantee  that 
the  equations  are  intrinsie  to  the  given  surfaees  of  interest.  We  will  further  verify  this  in  §2.4  to  help  the 
reader  with  the  intuition  behind  this  framework. 

®We  have  used  as  before  the  notation  A[x,^  =  A  x 
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2.2  Harmonic  Maps 

The  expressions  derived  in  the  previous  seetions  eome  from  the  theory  of  harmonic  maps,  e.g.,  [4,  7,  11, 
13,  15,  16,  20,  23,  35,  38,  39,  40].  In  general,  harmonic  maps  are  defined  as  maps  between  two  manifolds 
(Al,  g)  and  (Af,  h)  minimizing  the  energy 


E[u]  =  f  e[u]  (Wm  (13) 

JM 

where  in  local  coordinates  the  energy  density  e[u]  is  given  by 

A  1 

e[u]{x)  =  -  hij{u{x))  —  (14) 

We  have  used  Einstein’s  summation  here,  where  repeated  indices  indicate  summation  with  respect  to 
this  index,  together  with  the  usual  notation  for  tensors.^  When  both  the  domain  and  target  manifolds  are 
represented  explicitly,  the  classical  case,  the  Euler-Lagrange  equation  corresponding  to  this  energy  is  given 
by  (see  [38]) 


AW  +  r^(5)9“^^|^  =  0  (15) 

where  A_a4  is  the  Eaplace-Beltrami  operator  (reduced  to  the  regular  Eaplacian  for  the  case  of  flat  domain 
manifolds)  and  (ft)  stands  for  the  Christoffel  symbols  of  the  target  manifold  evaluated  at  u.  Note  that  the 
first  component,  the  Eaplace-Beltrami,  addresses  the  domain  manifold,  while  the  second  term  addresses  the 
target  manifold.  By  embedding  the  target  manifold,  we  are  changing  the  Christoffel  symbols  (expressing 
them  in  implicit  form,  see  below),^  while  the  work  in  [3]  changed  the  other  terms,  since  the  embedding  was 
done  to  the  domain  manifold,  see  §5. 

As  an  example,  let’s  see  what  happens  with  the  above  energy  for  the  Euclidean  case.  Since  both  metrics 
are  proportional  to  the  identity. 


e[u\{x) 


ij 


du^  \ 
dxj) 


2 


which  is  just  a  constant  multiplying  ||  Ja|l3^-  Therefore,  the  energy  defined  in  the  previous  case  is  just  a 
particular  case  of  harmonic  maps.  In  general,  this  energy  can  be  used  in  problems  such  as  color  image 
denoising  and  directions  denoising  [40,  41],  as  a  regularization  term  for  ill-possed  problems  defined  on 
general  surfaces  [17],  for  general  denoising  [37,  44],  for  models  of  liquid  crystals,  and  as  a  component  of  a 
system  for  surface  mapping  and  matching  [13,  31,  46]. 


2.2.1  An(other)  Informal  Calculation 

We  now  present  an  additional  computation  that  connects  in  a  deep  way  the  implicit  framework  with  harmonic 
maps.  We  consider  the  harmonic  energy  density  given  in  (14)  for  the  planar  domain  manifold  case  (p^  = 
5ij).  We  can  simplify  things  to  obtain 

e[u]{x)  =  ^  hij{u{x))  ^ 

V  V  p 

np“')o  =  A. 

®Or  alternatively,  the  second  fundamental  form  of  the  target  manifold. 
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We  know  that  Ily^  =  I  —  V-ip'V'ip'^  can  be  thought  of  as  the  inverse  of  the  target  manifold’s  metric 
tensor.  But  since  V'i/’  is  a  zero  eigenvalue  eigenvector  for  it  will  be  a  oo  eigenvalue  eigenvector  for 

n-J,.  Then,  we  can’t  use  the  identification  h  =  (hij)  -o-  in  the  above  expression  for  the  energy 

density.  However,  we  can  proceed  as  follows.  Take  e  >  1  and  define  fhe  mefric® 

h^=  (el- 

one  can  fhen  compufe  fhe  inverse  as  (if’s  an  elemenfary  formula,  see  for  example  [24]) 

h.  =  i(i  +  2^) 

The  energy  densify  can  be  rewriffen  as  (we  will  use  a  subindex  e) 

ee[u\{x)  =  ^  +  •VV’pj 

After  computing  fhe  variational  derivative  for  fhe  functional  [ft]  (x)  dx  we  obfain  fhaf  u  musf  satisfy 

Au  +  H.^[uxi,Uxi]  +  Au-V-ip^  V-ip  =  0 

By  multiplying  all  fhe  terms  in  the  above  equation  by  e  —  1  and  letting  e  — 1  we  find  fhaf  fhe  expression 
befween  brackefs  musf  vanish.  As  we  will  see  in  §2.4,  whaf’s  befween  brackefs  is  nofhing  buf  Av{x)  where 
oix)  =  'tp{u{x)).  So  is  a  harmonic  function  in  M.  If  is  also  evidenf  fhaf  o  satisfies  Dirichlef  boundary 
condifions  if  u  does,  and  since  we  are  frying  fo  map  fhings  from  Af  fo  J\f,  fhose  boundary  conditions  for  u 
musf  be  such  fhaf  dist{u{x),M)  =  0  for  x  €  M,  so  o\gj^  =  0.  Then  we  conclude  fhaf  o  musf  be  zero 
everywhere  in  Af . 

2.3  The  Mapping  Flow 

The  PDE  used  for  solving  fhe  harmonic  energy  is  given  by  ifs  corresponding  gradient  descent  flow.  This 
gradienf  descenf  is  given  by 


du^  i 

—  Au  +  H^(u) 

k=l 

where  fhe  inifial  dafum  uq  is  given  by  fhe  vecfor  field  we  wanf  fo  process,  fogefher  wifh  Neumann  boundary 
conditions: 


du  du  d'll^ 
dxk  ’  dxk  du^ 


{u) 


(16) 


f  u(a;,0)  =  uo(a^),  x^M 

\  iu^\dM  =  0.  ^  ^ 

The  use  of  Neumann  boundary  condifions  needs  fo  be  jusfified.  In  fhe  scalar  case,  one  has  fhe  evolufion 
problem 


It{x^t)  =  AI{x^t)  X  E  Ad,  f  >  0 
l(x,0)  =  Io(x),  xeM 
VI  ■  Ta\QM  =  0- 


® Since  e  >  1  all  the  eigenvalues  are  positive. 


(18) 
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We  observe  that  the  quantity  a{t)  =  I{x,  t)  remains  eonstant, 

a{t)  =  /  It{x,t)dMV=  /  AI{x,t)dMV=  /  V-{VI)dMV=  /  V7  •  n  =  0 
Jm  JM  JM  JdM 

thereby  imposing  the  boundary  eonditions. 

One  wonders  whieh  quantity  is  preserved  thru  time  by  the  flow  in  the  general  ease,  when  imposing  the 
boundary  eondition  (17).  We  illustrate  this  for  the  partieular  ease  of  AC  =  S^.  In  this  ease,  the  evolution 
equations  are  given  by  (see  also  §2.4  below) 

r  Xt  =  Aw  +  (||vxf +  ||vyf)w 

\  yt  =  Ay  +  diVATf +  ||vyf)y  ^ 

The  Neumann  boundary  eonditions  for  this  ease  are  written  as 

VA  •  n  =  Vy  •  n  =  0  inOM 

Transforming  to  polar  eoordinates  (p,  0)  one  finds  that  the  evolution  equations  (for  smooth  initial  data, 
and  at  least  for  some  time)  are  (see  also  [35]) 


with  boundary  eonditions 

Vd  •  n  =  0  in  dM. 

Again  one  finds  that  0{x,t)  dMV  is  constant. 

In  the  most  general  ease,  when  the  target  manifold  is  arbitrary,  one  might  guess  that  the  intrinsic 
barycenter^  of  the  map  is  preserved  through  time,  sinee  that’s  exaetly  what  the  partieular  eases  given  above 
show  us.  However,  to  the  best  of  our  knowledge,  there  is  not  sueh  a  result  in  the  literature  of  harmonie  maps, 
and  the  eonservation  of  the  baryeenter  is  only  obtained  when  eonstraints  are  added.  The  examples  diseussed 
above  still  motivate  the  use  of  Neumann  boundary  eonditions. 

2.4  Simple  Verifications 

We  now  illustrate  that  the  Euler-Lagrange  (12),  and  its  eorresponding  gradient  deseent  flow  (16),  are  the 
extension  for  implieit  targets  of  eommon  equations  derived  in  the  literature  for  explieitly  represented  mani¬ 
folds.  We  also  explieitly  show  that  the  flow  equation  guarantees,  as  expeeted  from  the  derivation  above,  that 
if  the  initial  datum  is  on  the  target  manifold,  it  will  remain  on  it.  We  also  express  the  seeond  fundamental 
form  of  a  manifold  that  is  implieitly  represented.  All  these  results  will  help  to  further  illustrate  the  approaeh 
and  verify  its  eorreetness. 

^The  intrinsic  baryeenter  G  of  the  map  ti  :  ^  W  is  defined  by  G  =  argminpes/^  /q  ttl*))  dx.  See  [10]  for  more 

details  on  the  baryeenter. 


Geodesics  as  Harmonic  Maps 

It  is  well  known,  see  [15,  16,  36],  that  are-length  parameterized  geodesies  on  the  manifold  Af  satisfy  the 
harmonie  maps  PDE.  If  we  assume  isotropie  and  homogeneous  metrie  over  Af,  we  end  up  having  that  (are- 
length  parameterized)  geodesies  must  satisfy 

7  +  [7, 7]  VV’(7)  =  0.  (21) 


Liquid  Crystals 


One  of  the  most  popular  examples  of  harmonie  maps  is  given  when  the  target  manifold  Af  is  an  hyper¬ 
sphere.  That  is,  the  map  is  onto  5*^“^.  In  this  ease,  the  embedding  (signed  distanee)  funetion  is  simply 
i’iy)  =  l|y||  -  y  e  From  this,  V-ipiy)  =  and  (H^(y)).^.  =  We  also  have  that 

since  ||5||  =  1.  In  addition,  u 


ply  obtained  taking  derivatives  with  respeet  to  Xk-  We  then  obtain  that 

,  2 


dxk  oxk  ^  J 


=  0,  faet  sim- 
2 

=  0,  and 


(fe”.)' 


from  (16)  is 


du  du 

dxif  ’  dxji 


Ylik  ( ~  ll*^u(^)  11^-  Therefore,  the  eorresponding  diffusion  equation 


du  ^  ^ 


ll- 


whieh  is  exaetly  the  well  known  gradient  deseent  flow  for  this  ease. 


Mapping  Restriction  onto  the  Zero  Level-Set 

We  now  explieitly  show  that  if  the  initial  datum  belongs  to  the  target  surfaee  given  by  the  zero  level-set  of 
then  the  solution  to  the  diffusion  flow  (16)  also  belongs  to  this  level-set.  This  further  shows  the  eorreetness 
of  our  approaeh. 

We  basieally  need  to  show  that  'tp{u{x,  t))  =  0,  Vx  G  AA,  Vt  >  0.  If  the  initial  datum  is  on  {-tp  =  0}, 
then  this  property  is  true  for  t  =  0.  Let’s  define  ^{x,  t)  =  'ip{u{x,  t)).  Then® 


dt 


„  , ,  du 

V^(«)  •  - 


d 

Au  •  V-ipiu)  +  ^  H^(w) 

k=l 


du  du 

dxk  ’  dxk 


sinee  is  a  distanee  funetion.  In  addition,  ^  =  V-tpiu)  ■  and  then 


d'^u 

dxf 


„  ,  du\  du  d'^u 


Adding  on  ?  =  1, ...,  d,  it  follows  that  ^  =  Ai^,  meaning  that  v  verifies  the  heat  flow.  In  addition  to  this, 
=  '^x  {'>P{u{x,t)))  ■  n  =  ■  n  =  (VV’(d))^Jan  =  {V'tp{u))'^0  =  0,  due  to  the  boundary 

eonditions  on  the  evolution  of  u. 

We  have  then  obtained  that  u  verifies  the  heat  flow  with  Neumann  boundary  eonditions  and  with  zero 
initial  data.  From  the  uniqueness  of  the  solution,  it  follows  that  u{x,  f)  =  0VxeAd,  Vf>0. 

®A11  the  calculations  that  follow  don’t  take  into  account  that  ip  might  fail  to  be  differentiable  at  some  points.  This  could  be 
addressed  by  a  regularization  argument. 
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Second  Fundamental  Form  for  Implicit  Surfaces 


If  we  compare  the  gradient  descent  flow  (and  Euler-Lagrange  equation)  we  have  obtained  with  the  classical 
one  from  harmonic  maps,  we  see  that  the  main  difference  is  that  the  Christoffel  symbols  of  the  target 
manifold  term  appearing  in  the  classical  formulation  has  been  replaced  by  a  new  term  that  includes  the 
Hessian  of  the  embedding  function.  We  obtained  this  by  first  embedding  the  target  manifold  and  then 
restricting  the  search  for  the  minimizing  map  to  the  class  of  maps  onto  the  zero  level-set  of  the  embedding 
function.  This  approach  can  be  followed  to  apply  this  framework  to  any  variational  problem.  We  now  show 
how  the  same  equation  can  be  obtained  by  simply  substituting  the  second  fundamental  form  of  the  explicit 
target  manifold  by  the  corresponding  expression  for  an  implicit  target  manifold.  This  will  illustrate  how  to 
apply  our  framework  to  general  PDE’s,  not  necessarily  gradient  descent  flow.  The  basic  idea  is  just  to  replace 
all  the  PDE  components  concerning  the  target  manifold  by  their  counterparts  for  implicit  representations. 

In  [28]  (page  150)  it  is  shown  that  the  scalar  second  fundamental  form  /i  at  a  point  p  of  an  hypersurface 
S  can  be  written  in  the  form 


h{p){V,W) 


^^{p)[v.w] 

||VV>P 


for  V,W  ^  TpS.  According  to  [28]  (page  139)  the  vectorial  second  fundamental  form  is  given  by 


I{p){V,W)=h{p){V,W)^^^ 

Prom  (15)  and  what  we  have  just  seen  it  is  obvious  that  the  implicit  version  of  the  harmonic  map  Euler- 
Eagrange  is  (12). 

As  stated  before,  and  following  the  formulas  above,  the  implicit  representation  of  the  target  surface 
permits  then  to  compute  the  second  fundamental  form  using  differences  on  Cartesian  grids,  without  the 
need  to  develop  new  numerical  techniques  on  polygonal  grids. 

Prom  the  result  just  presented,  in  order  to  transform  a  given  PDE  into  its  counterpart  when  the  target 
manifold  is  represented  in  implicit  form,  all  what  needs  to  be  done  is  to  re-write  all  the  characteristics  of  the 
PDE  concerning  this  target  manifold  in  implicit  form.  Por  completeness,  in  Appendix  1  we  present  basic 
facts  on  calculus  on  implicitly  represented  surfaces. 


2.5  Explicit  Derivation  of  the  Diffusion  Flow 

Here  we  first  proceed  in  a  naive  way  to  obtain  an  equivalent  formulation  of  the  gradient  descent  flow  that 
will  help  in  the  numerical  implementation.  We  assume  we  have  a  family  {ft (5,  t)}^of  mappings  from  H  to 
Af.  Por  each  t  we  define  the  harmonic  energy  of  a  member  of  the  family  as 

E{t)  =  ^  Jj\Ju{x,t)\\F  dx 

We  then  find  a  variation  of  the  family  such  that  E{f)  decreases.  To  accomplish  this  we  formally  differ¬ 
entiate  the  energy  with  respect  to  t.  A  simple  computation  yields 

E{t)  =  —  /  Uf  Au  dx 

Jn 

Now,  since  u{x,  t)  ^  J\f  y  x  €  H  and  V  f  of  smooth  existence,  one  must  have  ut{x,  t)  € 
appropriate  choice  for  ut  would  be 


ut  =  nT-(-,)^(A5) 


(22) 
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since  this  makes  E{t)  =  —  Jq  dx  <  0. 

The  projection  operator  in  (22),  as  we  already  know  (see  Appendix  1),  can  be  expressed  in  a  very  simple 
form  using  ip  (the  signed  distance  function  to  M), 

'^T^N'iv)  =  v  -  V  ■  Vipip)  Vip{p) 

Now,  it  should  happen  that  (22)  is  equivalent  to  (16).  We  show  this  in  §7. 


2.6  Remarks  on  the  Solutions  of  the  Diffusion  Flow 

The  well  posedness  of  the  diffusion  problem  with  Neumann  boundary  conditions  is  addressed  in  [22,  32], 
where  the  following  results  are  obtained,  here  included  for  completeness: 

Theorem  1  For  a  given  mapping  tto  :  At  — >■  AC  C  with  ^  =  0  on  dM.  and  for  every 

2  +  dim{A4)  <  p  <  +oo  there  exists  an  e  >  0  (depending  on  uq)  and  a  mapping  ft  :  At  — >■  AC  of  class 
L2(At  X  [0,  e],  ®  Moreover,  u  is  unique  and  C°°  except  at  the  corner  dM.  x  {0}. 

Theorem  2  Let  (M,g)  and  (AC, /i)  be  compact  Riemmanian  Manifolds  with  convex  boundary.  Let  u  : 
At  X  [0,a;)  Af  be  a  maximal  solution  of  the  diffusion  problem  with  initial  value  a  C°°  mapping  uq,^^ 

with  xo  =  ||e[uo]||Loo  >  0-  r  ^  Mbe  such  that  Ricm  >  <^^d  i?  >  0  such  that  all  sectional 

curvatures  of  Jf  are  not  greater  than  Then, 


1.  In  the  case  r  +  Rxo  >  0 


(a) 

(b) 


ut  > 


if  R>  0  then 
if  R  =  0  then  uj  =  +oo. 


w  >  i  log(l  +  when  r  ^0 
when  r  =  0 


2.  In  the  case  r  +  Rxo  <  0,  w  =  +oo. 


3  Maps  onto  Open  Surfaces 

So  far,  we  have  only  addressed  the  case  when  the  target  surface  is  closed  (zero  level-set).  In  this  section  we 
briefly  deal  with  open  surfaces.  We  show  that  when  the  function  is  evolving  according  to  the  flow  in  §2.3, 

the  set  C{t)  =  {u{x,t),  x  G  At}  remains  inside  the  initial  convex-hull  of  Co  =  {rto(a^),  x  G  At},  Vt  >  0. 
This  property  is  basically  a  consequence  of  the  maximum  principle.  Numerically,  this  might  of  course  be 
violated  due  to  numerical  errors,  and  we  will  later  discuss  how  to  correct  for  this  as  well. 

(At  X  [0,  e],  is  the  space  of  functions  /  :  At  — ^  such  that  for  every  i  =  1, . . .  ,n -\- 1,  VmR,  and  ^ 

are  all  in  L^(A(  x  [0,  ej). 

A  solution  it  :  At  x  [0,  oj)  — ^  A  of  the  diffusion  problem  is  maximal  if  it  cannot  be  extended  to  be  a  solution  on 
At  X  [0,  oj  -b  e)  for  any  e  >  0  or  if  oj  =  -boo. 

^^RicM  stands  for  the  Ricci  curvature  tensor  of  At. 
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3.1  Motivation:  The  Planar  Case 


Assume  that  the  target  manifold  J\f  is  flat,  for  example  (we  still  assume  that  the  domain  manifold  At  is 


dM 


=  0.  Let  S  be  a  eonvex  set  of  R 


k 


flat).  Let  u{x,t)  solve  ^  =  Art  for  x  €  Ad  and  f  >  0,  and 
with  smooth  boundary  (this  guarantees  that  the  distanee  funetion  is  also  smooth  almost  everywhere,  see  [36] 
for  a  formal  statement),  and  ^  the  signed  distanee  funetion  to  this  set  (positive  outside  and  negative  inside). 

Define  g{x,  t)  =  ^{u{x,  t)).  It  then  immediately  follows  that  ||  —  Ag  =  —  Yli=i  ^^(^5  Sinee  E 

is  eonvex,  so  it  is  Then,  the  Hessian  of  ^  is  positive  semi-definite,  meaning  that  ||  —  Ag  <  0.  Following 
the  sealar  maximum  prineiple,  max.[^£_;v(,t>o}  ■  If  x  €  Ad}  C  E, 

whieh  is  equivalent  to  0  >  ^{uo{x))  =  g{x,  0),  we  obtain  that  g{x,  t)  <  0,  and  u{x,  t)  G  E,  for  all  x  €  Ad 
yf  >  0. 


3.2  The  General  Case 

The  main  result  presented  below  is  from  [22].  We  quote  it  here  for  eompleteness. 

Theorem  3  Let  u{x,  t)  be  the  solution  of  (16)  at  time  t.  Let  us  assume  that  for  t  <T  this  solution  remains 
smooth.  Let  Lq  =  uo{Ll),  and  Iq  be  the  eonvex  hull  of  Lq.  Then  for  {x,  f)  G  x  [0,  T],  u{x,  t)  €  To- 


4  Maps  onto  Implicit  Submanifolds 

Here  we  present  a  modifieation  to  the  diffusion  flow  previously  presented  suited  to  diffuse  data  that  belongs 
to  a  eertain  submanifold  C  oi  N  =  {f)  =  0}.  We  speeify  the  submanifold  by  {■^  =  0}  fl  {$  =  0},  where 
we  seleet  $  :  IR^  — >■  i?  to  be  the  signed  intrinsic  (to  AC)  distanee  funetion  to  =  0},  satisfying  (see 
Appendix  1  for  the  notation) 


1  =  ||V^$||  =  VIIV^P  -  IVV”  V$|2  (23) 

In  addition  we  speeify  the  eondition 


$(p)  =  0  for  p  G  /C$ 

where 


/C$  =  {x  e  IR^  \  x  =  p-\-  aV'^(p),  with  p  ^  C,  a  ^  JR} 

is  the  cone  interseeting  {-ip  =  0}  at  C  and  direetor  rays  normal  also  to  {■^  =  0}. 

The  reason  for  speeifying  the  submanifold  this  way  is  that  we  eannot  proeeed  as  before,  simply  speeify- 
ing  the  submanifold  as  the  zero  level  set  of  it’s  Euelidean  distanee  funetion.  This  is  beeause  sueh  funetion 
would  be  singular  preeisely  onto  the  submanifold. 

As  we  show  in  Appendix  1,  the  Hessian  of  $,  intrinsic  to  AC  evaluated  at  the  point  p,  and  restrieted  to 
TpAC,  ean  be  written  in  the  form 


(p)  =  H,j(p)  -  A{p)  H^ip)  (24) 

where  A{p)  =  V$(p)  •'V'ip{p).  This  expression  will  be  used  below. 

^^Note  once  again  that  we  are  omitting  details  regarding  the  correct  handling  of  the  distance  function,  since  it  is  not  everywhere 
differentiable.  However,  by  a  regularization  argument,  the  same  conclusion  holds. 

^®The  proof  of  this  result  has  a  lot  of  interest  in  itself  since  it  can  be  carried  out  within  the  implicit  framework  introduced  in  this 
paper. 
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4.1  The  Minimization  of  the  Functional 

We  now  derive  the  Euler-Lagrange  eorresponding  to  this  additional  mapping  restrietion.  For  this,  we  use  a 
teehnique  slightly  different  that  the  one  in  §2.1. 

Let  us  assume  that  u  aehieves  a  minimum  of  the  energy  funetional  (1).  We  must  build  a  variation  of  u 
that  belongs  to  C,  the  interseetion  of  the  zero  level-sets  of  two  embedding  funetions  (and  not  just  of  as 
before).  It  is  elear  that  one  sueh  variation  would  be 


wx  =  lie  (ft  -f  An) 


We  are  interested  only  on  those  terms  of  ejwA]  that  do  not  vanish  after  the  g|“(*) 


A=o 


operation, 


namely  those  linear  in  A.  Therefore  we  only  preserve  those  terms  in  wx  whieh  are  eonstant  or  linear  in  A: 


We  write 


wa  ~  ft  +  AnT5c(i^) 


^Ti:c{v)  =  nr_{^=o}  {  ^  -  (n  •  V^$(u))  } 

=  V  —  {v  ■  V^$({t))  —  {v  ■  V'ip{u))  V'ip{u) 

where  —  A({t)  V'ip{u)  is  the  gradient  of  $  intrinsic  to  {-ip  =  0}. 

In  this  way  we  find  that  (up  to  a  first  order  in  A): 


N 


e[wA] 


e[u]  +  A  [vxi  -  Vxi  ■  V^$(u) 


(25) 


i—1 


Vxi  ■  yipiu)  V'tp{u)  —  V  •  — - - -V'tp{u)  —  V  •  V'tp{u)- 


dxi 


dxi 


Sinee$({t)  =  '4>{u)  =  0,  differentiating  with respeet  to  Xj  weobtainthat  =  V'tp{u)-Uxi  =  0, 

and  therefore 


-Uxi  =0 

The  expression  (25)  ean  be  simplified  fo  obfain 


N 

e[wx]  ^  e[u]  +  A  '^Uxi' 

i—1 


Vxi-V-  V^$(u) 


dxi 


V'ip{u) 


dV-ip{u) 

dxi 


Moreover,  sinee 


c)V^$(tt)  _  dA  A/-\xx  - 

- 

we  have 


dxi 


Ut,;  = 
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With  all  this  in  mind  we  find  that  (again,  up  to  first  order  in  A) 


N 

e[wx]  ^  e[u]  +  A  '^Uxi  ■  [vxi  ll!^[uxi,Uxi]  -v-V-ipiu)  ll^{u)[uxi,Uxi]] 

i—1 

Using  this  expression,  after  imposing  that  Jq  e[tyA]  dv  =  0  for  every  v,  we  find  fhaf  fhe  Euler- 

Lagrange  is 


Art  +  Eh, 


du  du 


dxk  dxk 


VV>(^)+ 


du  du 


dxk  dxk 


=  0. 


(26) 


an  expression  ufferly  prediefable. 


4,1,1  Simple  Verification 

As  for  fhe  ease  of  elosed  manifolds,  we  now  verify  fhaf  in  fael  fhe  gradienf  deseenf  eorresponding  fo  fhe 
Euler-Lagrange  (26)  keeps  ft  in  {■^  =  0}  fl  {$  =  0}.  We  jusl  need  fo  show  fhaf  i'{x,t)  =  'tp{u{x,t))  and 

IJ,{x,  t)  =  ^{u{x,  t))  are  always  zero.  The  idea  is  fhe  same  we  used  in  §2.4,  if  is  enough  fo  show  fhaf  bofh  v 
and  jjL  satisfy  fhe  heat  equation  wifh  adiabafie  boundary  conditions. 

1.  -ip 


We  have 


sinee  _L  Also 


and 


V-ip  •  Au  +  Eh* 

k 


du  du 
dxk  ’  dxk 


Au  =  V'tp  •  Au  +  Eh* 

k 


du  du 
dxk  ’  dxk 


ut  =  Au 


2.  $ 

We  have 


ht 


V$  •  Art  + 


du  du 


[dxk'  dxkj 


+  A 


du  du 


dxk  ’  dxk 


Erom  V$  •  =  ||V^$|P  =  1,  fhe  above  equation  eonfinues  as 
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V$-A5+ 

\  k 

V$ • Au +  (eh. 


du  du 


dxk  ’  dxk 
du  du  ’ 


+  A 

V  k 


dxk  ’  dxk 


Also 


and  then 


Ajj,  =  V$  •  Au  +  E 


Ht  =  Afi 


du  du 


dxk  ’  dxk 


du  du 


dxk  ’  dxk 


Finally,  it  is  easy  to  see  that  both  u  and  jj,  satisfy  Neumann  boundary  conditions.  Sinee  at  f  =  0  both 
funetions  are  zero,  we  must  have  that  they  are  identieally  zero. 


5  Implicit  Domain  Manifolds  and  p-Harmonic  Maps 

For  eompleteness,  we  present  now  the  formulas  eorresponding  to  the  ease  where  both  the  domain  and 
target  manifolds  are  represented  in  implieit  form  (with  the  implieitizing  funetions  being  the  eorresponding 
signed  distanee  ones).  Deriving  these  formulas  is  straightforward  using  the  framework  here  presented,  when 
eombined  with  the  work  in  [3].  We  also  show  the  eorresponding  flows  for  p-harmonie  maps. 


5.1  p-Harmonic  Maps 

We  still  assume  M.  to  be  planar.  The  energy  density  (2)  (but  no  the  dependenee  of  the  energy  on  its  density) 
is  redefined  as  follows.  For  every  p  E  [1,  +oo)  let 

ep[u]  =  -llJulljr 

A  simple  applieation  of  variational  ealeulus  leads  to  eonelude  that^^ 

Ut  =  p^~r  nv^(a)  (v  •  (^{ep[u\y~p  (27) 

Note  that  if  p  <  2  diffieulties  are  expeeted  to  arise,  see  [40]  and  the  referenees  therein. 


5.2  Generic  (Implicit)  Domain  Manifolds 

Let  M.  =  {  X  E  :  (j){x)  =  0},  where  (j){-)  is  the  signed  distanee  funetion  to  At,  then  the  diffusion  is 
given  by: 


k,r 


kr 


^^The  divergence  operator  convention  (for  a  matrix  A)  we  have  used  is  V  •  4  =  •  Ayi  |  . . 

the  i-th  column  of  A.  That  is,  we  apply  a  columnwise  divergence. 


VV’  (28) 

I V  •A,,.'),  where  Ayi  stands  for 
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The  whole  deduetion  rests  upon  the  redefinition  of  the  energy  (1)  and  its  density  (2).  Now  we  should  define 
fhe  energy  densify  fo  be 


eAu\  =  ^IIJ 


where  fhe  intrinsic  Jacobian  of  u  ean  be  wriffen  as  (see  Appendix  1  for  more  defails) 
The  new  definilion  for  fhe  energy  should  be:^^ 


E[u]  =  f  e(j,[u]  6{(f){x))  dx 

Comparing  fhis  wifh  (15),  we  ean  infer  fhe  implieif  form  of  fhe  Chrisloffel  symbols:^® 


(29) 


5.3  Generic  (Implicit)  Domain  Manifold  and  p-Harmonic  Maps 

Using  bofh  generalizafion  presenfed  above,  we  arrive  af  fhe  following  formula  wifh  a  bif  more  eompufafional 
efforl 


where 


=  p  n 


ut=p 


•Vip{u) 


(30) 


e</,,p[u]  —  ^1 


U  ''J- 


6  Diffusion  of  Tangent  and  Normal  Directions 


Throughouf  fhis  section  we  will  assume  dim{M.)  =  dim{JJ).  Assume  we  wanf  fo  diffuse  infrinsic  vectorial 
dafa  consfrained  fo  be  a  direcfion  (unif  norm)  and  fo  be  eifher  normal  or  fangenf  to  fhe  domain  manifold. 
We  can  fhen  minimize  fhe  functional  (29)  faking  a  variation  of  fhe  form  (assume  u  minimizes  fhe  energy 
functional  while  satisfying  bofh  ||{t||  =  1  and  If  (ft)  =  u) 

^  .  A  u  +  An(u) 

=  ||fl  +  An(j)|| 

where  v  :  A4  —>■  JR^  is  smoofh  and  If  is  eifher  or  ^NxM  (projecfion  onfo  fhe  fangenf  or  normal 

space  respectively).  Lef  w  =  n(u),  fhen  if  follows  easily  fhaf 


dE[ux] 


dt 


A=0 


—  /  {A^tt  +  2e^[u\  u}  ■  w  6{(f){x))dx 


Imposing 


A=o 


=  0  for  all  V  implies 


n  +  2  e^[u]  ft)  =  n  (A<^{t)  +  2  e^[u]  u  =  0 

^®We  have  already  taken  into  account  that  ||  Vi(>||  =  1. 

^®Of  course  pO  =  (=  Si)^)-  Then,  it  is  nice  to  observe  (although  formally  incorrect)  that  since  Ilv^  Vif>  =  0,  then  the 

metric  g  :  JR^  — ^  has  eigenvalue  +oo  in  the  direction  given  by  Vi(>  thus  prohibiting  intermingling  of  information  between 

adjacent  level  sets  of  (f). 
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Finally,  the  diffusion  flow  obtained  is 


—  =Ux{A^u{x,t))  +  2e^[u]{x,t)  u{x,t)  (31) 

Note  that  if  the  PDF  (31)  admits  a  smooth  solution  until  time  T,  and  if  (for  instance)  we  are  dealing 

with  tangent  directions  diffusion,  the  function  f{x,t)  =  V4>{x)  •  u{x,t)  satisfies  ft{x,t)  =  2e<^[{t]/(a;,  f). 
Therefore 

f{x,t)  = 

thus  verifying  that  if  V(f){x)  •  {t(a;,0)  =  0  then  V (f){x)  ■  u{x,t)  =  0  for  f  <  T.  We  also  want  to  check 
whether  ||{t(a;, f)||  =  1  V  (x^t).  Let 

Fci,[u]{t)  =  ^  [  llJalll-  dx 

then  Ffj,[u]{t)  =  —  ut  ■  A^u  6{(j){x))  dx.  Since  both  \\u\\  =  1  and  n({t)  =  u  (so  !!({<()  =  ut  since  If 
does  not  depend  on  t)  must  hold,  and  in  order  to  make  F^[u]{t)  non-positive  we  choose 


Ut  =  (32) 

where  nr„-p|a||=c}  =  I  “  MM  >  °- 

Note  that  the  above  evolution  indeed  forces  u{x,t)  to  satisfy  both  imposed  conditions.  Let  v  :  ]R^  — >■ 
]R^  be  such  that  n(n)  =  0  then  {v  •  u)t  =  v  •  ut  =  v  -  H  A^tt  =  =  0,  since 

the  projection  matrix  is  symmetric,  and  just  using  this  we  have  (|||f<|P)(  =  u  ■  Ut  =  u  ■  =  0 

trivially.  Finally,  using  ||{t||  =  1  and  carrying  out  some  computations  in  a  way  similar  to  §7  below,  one 
can  prove  that  (32)  reduces  to  (31). 

7  Numerical  Implementation 

We  now  discuss  the  numerical  implementation  of  the  flows  previously  introduced.  Since  the  target  manifold 
is  now  implicitly  represented,  we  can  basically  use  classical  numerical  techniques  on  Cartesian  grids.  Al¬ 
though  as  we  have  shown,  the  flows  guarantee  that  the  map  remains  on  the  target  (sub-)manifold,  numerical 
errors  can  move  it  away  from  it,  requiring  a  simple  projection  step. 

When  dealing  with  submanifolds,  although  the  evolution  equations  also  guarantee  that  the  solution  will 
remain  inside  the  convex  hull,  once  again  due  to  numerical  discretization  u  could  be  taken  outside  of  it  during 
the  evolution.  In  order  to  numerically  project  it  back,  we  need  to  have  a  distance  function  to  this  convex  hull 
defined  on  the  implicitly  defined  target  manifold.  In  [30]  we  have  shown  how  to  computationally  optimal 
compute  such  a  distance  function  on  implicitly  defined  manifolds,  and  this  is  the  technique  used  for  this 
projection  into  the  convex  hull. 

An  explicit  scheme  can  be  devised  to  implement  (28).  However  it  turns  out  that  it  is  more  convenient 
to  implement  a  mathematically  equivalent  evolution,  as  shown  in  [12].  More  specifically,  the  equivalent 
evolution  is 


—  =  Au  —  (Art  •  V'ip)'V'ip  (33) 

That  both  evolutions  are  equivalent  is  easy  to  see,  and  we  show  this  next. 

^^The  main  difference  is  that  now  one  must  take  into  account  the  Laplace-Beltrami  expressed  “implicitly,”  see  Appendix  1  for 
more  details  on  intrinsic  differential  operators  within  the  implicit  framework. 
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One  has  that  f{x,t)  =  'ip{u{x,t))  =  0  y{x,t)  e  0  x  iR+  U  {0}  for  u{-,-)  satisfying  (16).  Now, 
differentiating  /  with  respeet  to  Xi  we  obtain 


V^({t)  -Uxi  =0 

Differentiating  again  with  respeet  to  Xi, 

,  Kajj ]  +  V •  Ux^xi  ~  0 

Summing  for  all  i, 

'^'H.^[uxi,Uxi]  +  V-Ip-  Au  =  0 

i 

and  using  the  previous  expression  we  derive  (33)  from  (16). 

7.1  Numerical  Scheme 

All  the  eoding  was  done  using  Flujos  as  the  main  eore  (see  [19])  and  VTK  (see  [49])  for  visualization 
purposes.  All  the  examples  below  were  earried  based  in  equation  (30).  Its  numerieal  implementation  is 
straightforward  (at  least  when  p  =  2).  We  used  forward  time  diseretization  (explieit  seheme),  and  for  the 
spatial  diseretization,  we  used  the  following  well  known  reeipe.  To  spatially  diseretize 

ft{x,t)  =  V  •  {K{x)Vf{x,t))  (34) 

(K(a;)  is  a  symmetrie  positive  semi-definite  matrix),  we  eonsider  backward  approximation  of  the  divergenee 
and  a  forward  approximation  of  the  gradient.  Let’s  explain  how  this  applies  in  our  situation,  and  for  that  we 
assume  p  =  2  in  (30).  Then  the  equation  we  have  to  implement  is 


If  we  don’t  take  into  aeeount  the  outer  projeetion  matrix,  every  eoordinate  of  u  evolves  aeeording  to 

ul{x,t)  =  V  •  (nv<^(a,)Vu*(a;,f)) 

having  for  eaeh  eomponent  the  same  strueture  than  the  model  evolution  (34).  We  then  borrow  the  above 
diseretization  for  our  evolution.  If  we  eonsider  the  eoupling  among  different  tt*’s  imposed  by  the  projeetion 
matrix  we  see  that  we  still  preserve  numerieal  stability  sinee  this  matrix  is  positive  semidefinite  and 

has  speetral  radius  not  greater  than  1.^*  In  more  detail,  it  ean  be  shown  after  some  ealeulations  (see  [21,  42]) 
that  for  the  seheme  (p  now  denotes  a  position  over  the  grid) 

+  At  P(i^)(V-  •  (Q(p)V+i^) ) 

the  stability  eondition  is  of  the  form  (A  = 

,  .  .  _ 5(p) _ 

“  P,u  p(P(tt))  max{S'^(p),  Z?^(p)} 
or 

^®Note  that  ||tJ||^  >  Ilvvi  [tf;  1(1  =  ll^ll^  “  —  0  fo>'  tJ.  We  have  used  that  <p  isa  distance  function. 
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1 


A  < 


max, 


,P(PW) 


mm 

p 


!"{ 


S{p) 


maxis'^  {p),D^{p)} 


} 


where  p(P{p))  stands  for  the  speetral  radius  of  the  matrix  P(p),  S{p)  =  J2ij  (Qijip)  +  QijiP  ~  Axe,)), 
and  D{p)  =  {qij{p)  —  Qijip  —  Axe,)).  In  our  ease  we  may  admit  D{-)  to  be  small  eompared  with 
S{-)  (given  the  identifieation  Q  Hy,/))  when  Ax  is  small.  This  ean  be  easily  related  to  the  eurvatures  of 
{(^  =  0}  giving  a  eondition  on  the  sampling  of  the  distanee  funetion  (^)  representing  the  domain  manifold. 
This  eondition  mainly  means  that  we  require  a  fine  enough  sampling  as  to  guarantee  that  the  ehange  in  the 
normals  to  the  level  surfaees  of  (j)  is  small  between  adjaeent  grid  points.  This  eondition  is  obviated  when  the 
domain  manifold  is  planar.  So  the  stability  eondition  beeomes 


max„p(P(tt))  maXpS'(p) 

Sinee  by  Cauehy-Sehwartz’s  inequality  (and  the  aforementioned  assumption  on  the  ehange  of  between 
adjaeent  grid  points)  2  d  (in  praetise)  upper-bounds  S{p),  remembering  the  faet  that  p(P(p))  <  1,  we  arrive 
at  A  <  .  Note  that  if  a  more  eareful  implementation  is  desired,  good  ehoiees  are  ADI  or  AOS  sehemes, 

see  [50] 

All  derivatives  in  and  were  approximated  by  eentral  differenees.  An  interpolation 

seheme  had  to  be  used  sinee  the  evaluations  of  in  the  above  equation  are  at  positions  given  by 

u{x,  t),  positions  not  neeessarily  on  the  underlying  grid.  We  used  linear  interpolation  for  this  purpose. 

Note  that  as  done  in  [3],  when  the  domain  manifold  is  also  implieitly  represented,  the  values  of  the  map 
on  it  are  periodieally  extended  to  its  surrounding  offset  due  to  stability  eonsiderations.  Also,  as  explained 
before,  due  to  numerieal  diseretization,  the  diseretely  eomputed  solution  map  ean  be  taken  out  of  the  target 
manifold  during  the  evolution.  In  this  paper,  we  simply  projeet  it  baek  at  every  iteration.  We  have  seen 
that  this  projeetion  is  a  trivial  step  due  to  the  faet  that  the  embedding  is  a  distanee  funetion.  It  is  quite 
straightforward  to  show  that  the  results  reported  in  [1]  ean  be  extended  for  our  equations  as  well,  at  least  for 
eonvex  hyper-surfaees. 


7.2  Numerical  Examples 

In  all  the  examples  below,  the  domain  manifold  M.  is  either  the  Euelidean  spaee  or  an  implieit  torus. 
The  target  manifold  Af  is  an  implieit  surfaee  in  JR^,  that  is,  the  zero  level-set  of  'll;  :  ]R^  — iR,  being  a 
signed  distanee  funetion  (this  is  of  eourse  also  the  ease  when  the  surfaee  is  a  sphere,  ip  being  as  in  §2.4). 

In  order  to  present  interesting  examples  we  eonstruet  texture  maps,  add  noise  to  them,  and  then  diffuse 
them  using  our  framework.  Let  S  be  the  surfaee  onto  whieh  we  want  to  map  a  given  (planar)  image  defined 
in  a  subset  D  C  IE?.  Then  the  texture  pap  is  a  map  T  :  S  ^  D.  Onee  the  map  is  known,  we  inverted  it  to 
find  a  map  'Uq  :  D  ^  S.  Then,  we  builf  up  fhe  noisy  map  u  ■.  D  ^  S  defined  by 

‘U,{x)  =  {‘u,q{x)  -h  n(a;)) 

where  n  :  D  — 5  is  random  map  wifh  small  preseribed  power  a.  We  fhen  feed  fhe  evolution  (16)  wifh  ft  as 
initial  eondifion,  and  Neumann  boundary  eondifions.  Afler  a  eerfain  number  of  sfeps  we  stop  fhe  evolution, 
inverf  fhe  resulfing  map,  and  use  if  as  a  fexfure  map  to  painf  fhe  surfaee  wifh  a  eerfain  fexfure. 

As  a  means  of  finding  a  suifable  T  we  have  exfended  fhe  work  in  [47]  (a  multidimensional  sealing 
approaeh),  eombined  wifh  fhe  feehnique  developed  in  [30]  for  eompufing  disfanees  on  implieif  surfaees. 

^®Note  that  we  are  not  proposing  this  as  a  complete  texture  mapping  alternative,  it  is  just  to  provide  an  illustrative  example. 
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Figure  1:  Diffusion  of  a  noisy  texture  map  (left)  onto  an  implicit  sphere  (right).  (This  is  a  color  figure.) 


In  all  the  steps  just  described  there  are  some  minor  implementation  details,  mainly  regarding  interpolation 
tasks,  that  we  omit  for  the  sake  of  clarity. 

In  Figures  1,  2  and  3  we  then  denoise  vectors  from  the  plane  ]R^  to  a  3D  surface  defined  as  the  zero 
level-set  of  t/)  :  IB?  M  and  map  a  texture  image  to  the  surface  using  the  obtained  map.  Note  that  the  map 
is  the  one  being  processed,  not  the  image  itself. 

We  also  show  an  example  of  diffusion  of  random  maps  from  an  implicit  torus  to  the  implicit  bunny 
model,  see  Figure  4.  As  expected  from  the  theory,  when  evolving  this  set  with  the  harmonic  flow,  the  set 
converges  to  a  unique  point. 

8  Conclusions 

In  this  paper  we  have  shown  how  to  implement  variational  problems  and  partial  differential  equations  onto 
general  target  surfaces.  We  have  also  addressed  the  case  of  open  target  surfaces  and  sub-manifolds.  The 
key  concept  is  to  represent  the  target  (sub-)manifolds  in  implicit  form,  and  then  implement  the  equations 
in  the  corresponding  embedding  space.  This  framework  completes  the  work  with  general  domain  mani¬ 
folds  reported  in  [3],  thereby  providing  a  complete  solution  to  the  computation  of  maps  between  generic 
manifolds. 
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Figure  2:  Diffusion  of  a  noisy  texture  map  onto  an  implicit  teapot.  We  show  two  different  views  (noisy  on  the  top 
and  regularized  on  the  bottom).  (This  is  a  color  figure.) 
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Figure  3:  Diffusion  of  a  texture  map  for  an  implicit  teapot  (noisy  on  he  top  and  regularized  on  the  bottom).  A  chess 
board  texture  is  mapped.  (This  is  a  color  figure.) 
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Figure  4:  Diffusion  of  a  random  map  from  an  implicit  torus  to  the  implicit  bunny.  In  blue  are  marked  those  points 
of  the  bunny’s  surface  pointed  by  the  map  at  every  instant.  Different  figures  correspond  to  increasing  instances  of  the 
evolution,  from  top  to  bottom  and  left  to  tight.  We  show  the  map  at  17  of  100  iterations  performed  to  the  initial  map 
with  a  time  step  of  .01.  We  used  the  2-harmonic  heat  flow  with  adiabatic  conditions.  (This  is  a  color  figure.) 
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Appendix  1:  Implicit  Calculus 


We  now  present  basie  faets  about  differential  ealeulus  on  implieitly  represented  surfaees.  For  more  infor¬ 
mation  see  for  example  [2,  8,  29]. 

We  have  a  smooth  sealar  funetion  /  :  ]R^  M,  and  a  smooth  veetor  field  A  :  ]R‘^  {d  and  D  are 

not  neeessarily  equal).  The  manifold  onto  whieh  the  ealeulus  is  to  be  done  is  represented  as  <S  =  {■^  =  0}, 
for  -^(•)  the  signed  distanee  funetion  to  S. 

All  the  ideas  of  differentiation  ean  be  obtained  from  simple  eonsiderations  related  to  the  restrietion  of  the 
funetion  to  a  geodesie  eurve  living  in  the  manifold.  We  eonsider  an  are-length  parameterized  geodesie  eurve 
7  :  [— e,  e]  — >■  5  sueh  that  7(0)  =  p  is  a  given  point  of  S.  We  denote  F{t)  =  /(7(f))  and  A(f)  =  X{'y{t)). 

Implicit  gradient 

We  differentiate  onee  F{t)  to  obtain  i^(0)  =  V  f{p)  ■  7(0).  Sinee  7(0)  G  TpS  (the  tangent  plane),  we  find 
the  implied  gradient  of  /  at  p  to  be  V5/(p)  =  V f{p)  —  V f{p)  •  n{p)  n{p),  where  n{p)  stands  for  the 
normal  to  the  manifold  atp.  Sinee  we  ean  also  write  n{p)  =  V'tp{p),  we  obtain 

V5/(p)  =  V/(p)  -  (V/(p)  •  v^(p))  v^(p) 

We  often  use  the  alternative  notation  V-ipf  sinee  the  definition  ean  be  applied  to  any  level  set  of  ■)/).  Note 
that  we  ean  write  V^/  =  IIv^V/  where 

=  I  - 


Implicit  Hessian 

If  we  eompute  the  seeond  derivative  of  F  we  find  fhaf  i^(0)  =  Vf{p)  ■  7(0)  -f  Hy[7(0), 7(0)].  Now, 
we  know  fhaf  an  are-lengfh  parameferized  geodesie  eurve  of  S  musf  safisfy  fhe  harmonic  maps  differenfial 
equafion 


7-f  H^(7)[7,7]  V^(7)  =  0 

We  fhen  find  fhaf  F{0)  =  (Hj(p)  —  V f{p)  ■  V'ip{p)  H^(p))[7,7].  Again  we  have  fhaf  7  G  TpS,  and  we 
find  fhe  implieif  Hessian  of  /  af  p  fo  be 


Hf(p)  =  n^hyn^ 

where 


hy  =  Hy(p)-V/(p)-VV>(p)  H^(p) 

We  will  frequenlly  use  fhe  alfernafive  nofafion  Hj  (p). 

Implicit  Laplacian 

From  fhe  previous  eompufafion  if’s  an  easy  exereise  fo  eompufe  fhe  implieif  Laplaeian  or  Laplaee-Belframi 
of  /  sinee  by  definifion  A5/  =  frace{Hj}. 

For  any  pair  of  symmefrie  mafriees  A  and  B  one  has  fhaf  frace{ ABA}  =  Yj  Yk  O'ijdikbjk  and 
frace{ AB}  =  Yi  have  fhaf  n^BII^  =  B  -|-  V-ipVip'^'BV'ipV'ip'^  —  VipV'ip'^'B  — 

'SS/ il)V ijF .  We  fhen  obfain 
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tracejn^Bn^}  =  trace{B,}  +  '^'^'^'ipxi'tpxj'tpxi'tpx^bjk 

i  j  k 

~  2  'tpxji’xjbjj 

1  3 

Recalling  that  '^(•)  is  a  distance  function,  so  that  it  satisfies  HV-^H  =  1,  we  find 

fracejn^BII^}  =  frace{B}  —  EE  'tpxi'fxjbij 

i  3 

=  tracefB}  —  B[V'^,  V-^] 

We  conclude  the  reasoning  by  taking  B  =  h j : 

trace{Hj}  =  fracejhj}  —  hj-[V'^,  V'^] 

=  frace{hj-}  —  Vt/)] 

since  V'i/)]  =  0.  Since  tracefUf}  =  Af  —  (V/  •  V'tp)A'tp,  we  find  that 

Asf  =  Af-  (V/  •  -  Hy[VV',  Vfj] 

It’s  interesting  to  observe  how  the  expression  just  found  for  Asf  coincides  with  the  one  obtained  by 
minimizing  the  intrinsic  Dirichlet  integral, 

D{f)  =  l[  \\VsffS{f;)dv 

2  Jm.d 

as  is  done  in  [3].  The  authors  showed  that  a  smooth  function  /  extremizing  D{f)  must  satisfy 

V  •  (V/ -  (V/ •  VV>)  Vf3)  =  0 

We  should  verify  that  this  definition  coincides  with  ours.  This  is  accomplished  as  follows: 

V  •  (V/  -  (V/  •  V^)  Vi>)  =  Af-  (V/  •  V^)A^  -  V(V/  •  V^)  • 

=  A/-(V/-VV>)AV>-Hy[VV>,VV>]  -H^[V/,VV>] 

=  A/- (V/- V^)A^-Hj[V^,V^] 

=  Asf  (according  to  our  definition), 

since  H^[VV>,*]  =0. 

^°As  one  expects  since  this  is  the  definition  of  harmonic  functions. 
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Vector  Calculus 


•  Implicit  Jacobian:  With  the  ideas  developed  before,  we  easily  find  (differentiating  A(f))  that 

•  Implicit  Divergence:  Using  the  expression  for  the  intrinsie  Jaeobian  we  write 

V5  •  A  =  trace 


and 

V5  •  A  =  V  -  A- J^[V^,V^] 

It  is  useful  to  observe  that  V5  •  A  =  V  •  A  when  X{x)  E  =  0} 
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